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Abstract 

Modelling the Calvin cycle of photosynthesis leads to various systems 
of ordinary differential equations and reaction-diffusion equations. They 
differ by the choice of chemical substances included in the model, the 
choices of stoichiometric coefficients and chemical kinetics and whether 
or not diffusion is taken into account. This paper studies the long-time 
behaviour of solutions of several of these systems, concentrating on the 
ODE case. In some examples it is shown that there exist two positive 
stationary solutions. In several cases it is shown that there exist solutions 
where the concentrations of all substrates tend to zero at late times and 
others (runaway solutions) where the concentrations of all substrates in- 
crease without limit. In another case, where the concentration of ATP is 
explicitly included, runaway solutions are ruled out. 

1 Introduction 

Photosynthesis is a process which is of great importance for many reasons. It 
is the ultimate source of the food we eat, the oxygen we breath and many 
fuels (fossil fuels and biofuels). For this reason it is clear that it would be 
valuable to have a better theoretical understanding of this process and one way 
of approaching this task is to use mathematical models. The aim of this paper 
is to analyse dynamical properties of some of these models. 

Photosynthesis can be split into two main parts. In the first part, called the 
light reactions, energy is captured from light and the small molecules ATP and 
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NADPH are produced. These provide sources of energy and reducing power, 
respectively, for the second part, the dark reactions. The name of the latter 
comes from the fact that they can take place in the dark. Molecular oxygen is 
produced during the first part. In the second part carbon dioxide is used to make 
carbohydrates. For this reason this part is also known as carbon fixation. This 
paper is exclusively concerned with models for the second part of photosynthesis. 

The models which are relevant describe reactions between different chemical 
substances and also, in some cases, diffusion of these chemicals. The resulting 
mathematical model is a system of ordinary differential equations (ODE) if 
diffusion is not included and a system of reaction-diffusion equations if it is. 
The models studied in what follows are either taken from the papers [S] or [52] 
or are closely related to the models in those papers. In all these cases the network 
of reactions modelled contains a cycle and due to the fundamental contributions 
of Melvin Calvin to identifying the reactions concerned this is often referred to 
as the Calvin cycle (see for instance [T]). 

In building a model it is necessary to decide which substances are to be 
included. The basic unknowns are the concentrations of these substances. It 
is also necessary to decide how the chemical reactions are to be modelled. In 
most of this paper diffusion is ignored and only a few remarks are made on what 
happens when it is included. In the absence of diffusion the equations are of the 
general form 

X, = /,;(x). (1.1) 

Here Xi are the concentrations, which are functions of time, and the dot denotes 
the time derivative. The solutions of relevance for the applications are those 
for which all Xi are positive. In other words the point with coordinates Xi{t) 
is always in the positive orthant S of R". The mapping / with components fi 
represents the interaction between the different substances during the reactions. 
It is of the form Nv{x) where is a matrix called the stoichiometric matrix and 
the components Va of v are the rates of the different reactions. The Va describe 
what is called the kinetics. The function fi is of the form — Xif~ for two non- 
negative functions and /~ as a result of the form of the dependence of the 
reaction rates on the concentrations of the substances going into the reactions. 
The cosets of the range of the stoichiometric matrix are called stoichiometric 
compatibility classes and are invariant under the flow of the dynamical system. 
In all the systems considered in what follows the function v is and it can 
be shown that if the concentrations Xi are positive at some time they remain 
positive as long as the solution exists. This can be proved as in the special case 
covered by Lemma 1 of |19) . Thus S is invariant under the evolution and it 
follows by continuity that its closure 5' is also invariant. 

A common choice of kinetics is mass action kinetics where if p molecules of 
the substance with concentration Xi take part in a reaction the reaction rate 
has a factor proportional to x^. This corresponds to the idea that the rate of 
reaction is proportional to the probability of the relevant molecules meeting. 
For instance in the simple reaction A + 2B C the reaction rate is of the form 
kxAx\ where k is the reaction constant. For more details on building systems 
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of ODE describing reaction networks and mass action kinetics in particular 
see [7]. Another common choice, particularly in the description of biological 
systems, is Michaelis-Menten kinetics. This is adapted to describing reactions 
which are dependent on a catalyst and the reactions in biological systems are 
usually catalysed by enzymes. If the enzymes are explicitly included in the 
description a type of kinetics is obtained which is referred to in [5] as Michaelis- 
Menten represented in terms of mass action (MM- MA) . Michaelis-Menten (MM) 
kinetics is obtained from MM-MA kinetics by a limiting process (quasistationary 
approximation) . 

The structure of the paper is as follows. In Section [2] the dynamics of mod- 
els with mass action kinetics is considered for two different choices of the sto- 
ichiometric coefficients. In particular it is shown that for certain values of the 
reaction constants there is exactly one positive steady state and that it is un- 
stable. This raises the question of the final fate of general solutions. It turns 
out that for suitable choices of the reaction constants there is an open set of 
initial data for which all concentrations tend to zero at late times and an open 
set of initial data for which all concentrations tend to infinity at late times. The 
main results are collected in Theorem 1. The second statement is somewhat 
technical to prove for the choice of stoichiometric coefficients used in [3] and 
the proof is the subject of Section [3] Section H] is concerned with the models 
where the kinetics is Michaelis-Menten represented in terms of mass action. It is 
shown that there are solutions which tend to infinity at late times for both the 
choices of stoichiometric coefficients made in [9] and in [22]. For the first case 
it is proved that there can exist more than one positive stationary solution in a 
stoichiometric compatibility class and there is some discussion of what happens 
in the second case. The models with Michaelis-Menten kinetics are studied in 
Section [S] and it is shown that that there are solutions which tend to infinity 
at late times for that model too. It is shown that the stationary solutions are 
essentially the same for the MM-MA and MM models. A model in which the 
concentration of ATP is a dynamical variable is discussed in Section [51 This 
leads to a system of ODE for which, in contrast to the models discussed up 
to this point, all solutions are bounded in the future. In all this the aim is 
to treat values of the reaction constants which are as general as possible. In 
Section [7] some conclusions are drawn. Appendix A gives an introduction to 
Michaelis-Menten theory. Appendix B collects some technical results required 
for the proofs in the main text. 

2 Mass action kinetics 

This section is mainly concerned with the dynamical system (6) of |^. There 
are five variables a;RuBP, a^pcA, a;DPGA, a;GAP and xrusp which are the con- 
centrations of the substances abbreviated by the subscripts. They are ribulose- 
1,5-bisphosphate (RuBP), 3-phosphoglycerate (PGA), 1,3-diphosphoglycerate 
(DPGA), glyceraldehyde-3-phosphate (GAP) and ribulose-5-phosphate (Ru5P). 
This system has mass action kinetics and is called MA in what follows. It is 
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given by 



BP 



dt 

dxpGA _ 

dt 

dxuPGA 

dt 



= hXRubP - fcl^RuBP, (2.1) 

= 2A:i2:ruBP - fc2a;pGA - hxpGA, (2.2) 

= ^22^PGA - ^3a;DPGA, (2.3) 



^^dt^ = fcaa^DPGA - 5kiX^j^p - fc7a:GAP, (2.4) 

rfXRu5P 



dt 



fcsa^RuSP + 3/c42^GAP- (2-5) 



The hi are the reaction constants and they are ah positive. The ahcrnative nota- 
tion where (xruBP, a;pGA, xbpga, a::GAP, a^RuSp) is replaced by (xi, X2, X3, X4, X5) 
is also used. The state space of interest for the applications is the positive or- 
thant S. Sometimes it is also useful to consider the dynamics on S. The origin 
is a stationary solution. The linearization of the system at the origin has eigen- 
values (— fci,— ^2 — kg, —k^, —kf, —k^). Thus the origin is a hyperbolic sink. 
Consider a solution which starts at a point of the boundary of S other than the 
origin. Let N be the set of indices i for which the concentration Xi vanishes. 
Both N and its complement are non-empty. Hence there exists i ^ N ioi which 
j & N for j = i + 1 mod 5. It follows that Xj > and so the extension of the 
solution towards the past must lie in the complement of S. This implies that 
there exists no solution other than the zero solution which stays in the boundary 
of S for a finite time. This can be used to show that if x is a solution which 
starts in 5* then its oj-limit set contains no point of the boundary of S other 
than the origin. For suppose that x* is a point of the w-limit set of x{t) which 
belongs to the boundary of S and is not the origin. Then there is a solution y 
which passes through x* and lies entirely in the w-limit set of x and hence in S. 
On the other hand it has just been shown that this cannot happen. Thus any 
oj-limit point of a solution starting in S must either be the origin or a point of 
5*. 

Taking a suitable linear combination of (|2.4p and (|2.5|) eliminates the non- 
linear terms. 

d(3xGAP + 5XRu5p) „, „, ^, 

— = 3fc3XDPGA - 3fc7XGAP ~ OfcsXRuSP- (2.6) 

Let X be the maximum of the quantities xruBP, xpga, xdpga, xgap and 
3xGAP +5xRu5P . Then any solution of the system satisfies the integral inequality 

X{t) < X{to) + C [ X{s)ds (2.7) 

J to 

where C is the maximum of 2fci, ^2, 3^3 and ^k^. Thus, by Gronwall's inequal- 
ity, none of the variables can blow up in finite time. Together with the fact 
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that S is invariant this shows that the solution exists globally in the future. 
Summing up what has been proved so far gives: 

Proposition 1 A solution of (|2.ip - (|2.5p with positive initial data exists globally 
to the future, remains positive and has no w-limit points on the boundary of S 
except possibly the origin. 

It is shown in j5] that if k2 > 5fc@ there is a unique stationary solution of 
the system in S and the equilbrium concentrations are calculated explicitly in 
terms of the reaction constants. For this solution 



a^RuBP = fci ^ 



3k A 



2k2 5 



k2+ke 3 



(2.8 



The other equihbrium concentrations can be expressed as xru5p — IJ-xruBP, 

^PGA = fc;+fe^a;RuBP, XDPGA = Mfc^+fcJ) ^^RuBP , XgAP = (^gfeja^RuBpJ • An 

additional relation which can be derived for the stationary solution is that 
■''GAP = kj(k2-5ks) • "^^^ linearization of the system about the stationary so- 
lution has the characteristic polynomial 

(A + fci)(A + fc2 + ke){\ + fc3)(A + 25/c4a:QAP + k7){\ + k^) - 'iQkik2ky,kikzXQp^^ . 

(2.9) 

The constant term is equal to 

kik'ik^[{k2 + ke)k7 + 5(-fc2 + 5kQ)kiXQp^p] = -4fci(A;2 + kQ)k^k^k-j. (2.10) 

Because of the signs of the coefhcients this polynomial has exactly one positive 
root. This means in particular that the stationary solution is unstable. For 
^2 < Sfce there is no stationary solution in S. If k2 — S/cg is allowed to tend 
to zero while each of the reaction constants tends to a non-zero value then the 
stationary point tends to infinity. 
Define a function 

13 3 , , 

Li = XRuBP + 2^PGA + gaiDPGA + g^GAP + a;Ru5P- (2-11) 

Then 

^^2 1 5 / ^^^^ ~ -Kya^GAP- (2.12) 

If k2 < 5^6 then L\ is a Lyapunov function for the system (|2.ip - (l2.5p . This 
recovers the fact that for this parameter range there are no stationary solutions. 
It fact, when combined with Proposition 1, it shows that when that inequality 
holds all solutions converge to the origin as i — >■ oo - all concentrations go to zero 
at late times. Thus strong control of the late-time asymptotics of all solutions 
has been obtained in this case. 

It remains to consider the case k2 > Sfcg where Li does not seem to give an 
interesting conclusion. A useful generalization of Li is given by 

L2 = a;RuBP + 2'^PGA + aa;DPGA + aa:^GAP + a;R,u5P (2-13) 
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for a > 0. It satisfies 

= [h - (2a - 1) k2] xpGA - [akj - (3 - 5a)kiXQj<^p]xGAP ■ (2.14) 

This means that if (2a — l)k2 < fcg the function L2 is monotone decreasing 
along any solution for which (3 — 5a)fc4XQAp < akj. For a = f the function 
L2 coincides with Li. Another interesting choice is a = i. In that case L2 is 
monotone decreasing on the region defined by the inequality k4XQji^p < kj. It 
follows that if a solution initially satisfies 



1/ 

a;RuBP + -^[XPGA + a;DPGA + xqap, 



1 fkj 



(2.15) 



then it tends to zero as t — > 00. 

It has now been shown that in the case k2 > S/ce there is an open set of 
initial data for which the corresponding solutions tend to the origin as i 00. 
The argument just given also provides some information about the basin of 
attraction of the origin, which is more than could be concluded from the fact 
that the origin is a hyperbolic sink. There is also an open set of initial data 
for which all concentrations Xi tend to infinity as i — >■ c». The proof of this 
statement is given in Section [31 

In [22] stoichiometric coefficients are considered which are slightly different 
from those in [9] . While [22] uses Michaelis-Menten kinetics it is possible to take 
mass action kinetics with the stoichiometric coefficients of [2^ . This leads to a 
system which is called MAZ in what follows. It differs from the system MA only 
by the facts that the terms — 5fc4a;Q^p and 3fc4a;Q^p are replaced by —kiXQAP 
and |fc4a;GAP respectively. In terms of the reaction network, the system MA 
arises from the system MAZ by multiplying the stoichiometric coefficients in 
one of the reactions by a constant factor so that they become integers. For the 
system MAZ the set S is positively invariant. The right hand side of the system 
is due to the fact that all the stoichiometric coefficients on the left hand sides 
of reactions are integers. The arguments used to prove that solutions starting 
in S have no cj-limit points on the boundary of S other than the origin for the 
system MA generalize easily to give the same statements for the system MAZ. 
Since the latter system is linear all solutions exist globally. The function Li of 
(|2.1ip is a Lyapunov function for the system MAZ when k2 < Sfcg since it also 
satisfies the equation (|2.12p in this case. It follows that all the statements about 
the system MA derived using Li also hold for the system MAZ. In the latter 
case 



^ = -2[^« 



{2a - 1) k2] xpGA 



ak^ 



l^a 



XGAP- 



(2.16) 



Thus if {2a — l)k2 > k^ and (3 — 5a)fc4 > 5aA;7 and at least one of these two 
inequalities is strict then —L2 is a Lyapunov function. The equations for sta- 
tionary solutions may be analysed in this case in a similar way to what was 
done for the system MA. Four of the relations obtained are xrusp = ^xp^^^p, 
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2;pGA = fell^XRuBP, a;DPGA = k.fkl+ke) ^^^^^^ ' ^GAP = Ija^RuBP. Substitut- 
ing these relations in the remaining equation gives: 

5{k2 + ke){ki + kj) = 6k2ki. (2.17) 

Thus there are stationary solutions only when this relation is satisfied and when 
it is satisfied there is a whole one-dimensional subspace of them. This situation 
is not surprising since the equations are linear in this case. It is natural to 
examine the eigenvalues of the matrix on the right hand side of the equation. 
The characteristic polynomial is 

6 

(fcl + A)(/C3 + A)(fc5 + A)(fc2 +k6 + X){k4 + kT + X) - -kik2k3kik5. (2.18) 

5 

As in the case of the system MA all the coefficients in this polynomial are 
positive except possibly for the constant term, which is 



kiksk^ 



6 

(fc2 + k(i){k4 + kr) - -^2^4 
5 



(2.19) 



When the expression in brackets is negative there is a positive real eigenvalue 
and there exists at least one solution which tends to infinity as < — > cx) since the 
origin has a non-trivial linear unstable manifold. If 

6 

{k2+kG){ki + k7) < -k2ki (2.20) 
5 

then a can be chosen so that —L2 is a Lyapunov function. Hence in that case 
the origin does not belong to the w-limit point of any solution. Since other 
w-limit points on the boundary of the positive orthant have already been ruled 
out it follows that all solutions tend to infinity. The results for the systems MA 
and MAZ are now collected as a theorem. 

Theorem 1 If ^2 < Sfcg then all solutions of MA and MAZ tend to the origin 
as t ^ 00. If k2 > Sfcg there is a non-empty open set of initial data for MA for 
which the corresponding solutions tend to the origin as t — >■ cxd and a non-empty 
open set of initial data for which the corresponding solutions tend to infinity as 
t — > 00. If ^2 > Sfcg then there is at least one solution of MAZ which tends to 
infinity as i — > 00. If /c4(fc2 — Sfcg) > 5(fc2 -I- ^6)^7 then all solutions of MAZ tend 
to infinity as t — >■ 00. 



3 Solutions which tend to infinity 

In this section it is shown that if k2 > Sfcg there exist solutions of (I2.ip - (j2.5p for 
which all Xi tend to infinity as t — ^ 00. In fact there is a non-empty open set of 
initial data for which the corresponding solutions behave in this way. 
Theorem 2 If ^2 > Sfcg there is a non-empty open set of initial data for the 
system MA for which the corresponding solutions tend to infinity as t — >■ 00 and 
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have the asymptotics 



= Ae"* + ... 
X2{t)^ A{V2/Vi)e^' + .. 
x:i{t) = A{V^/V^)e'^' + .. 
Xi{t) = [A{V:i/V,)i,]^e^ 
X5(i)=^(14M)e"* + .. 



(3.1) 
(3.2) 
(3.3) 
(3.4) 
(3.5) 



Here a, and the Vi are fixed positive constants and ^ is a constant depending 
on the solution. 

Proof For the proof it is useful to introduce the quantity — —■ Its evolution 
equation is given by 

5(a;30^ 



dt 



♦ / 1 , , k2X2 

-ks - kj - - — 

5 0x3 



(3.6) 



Let X be the vector with components Xi for i = 1,2,3,5. Then four of the 
evolution equations can be rewritten as 



dx 

— = Mx- 
dt 



R. 



(3.7) 



Here the only non-zero component of the vector R is the last one and it is equal 

to 3/04x3(^-6), 6 = ^ 



M 



and 








-fci 










2fci 


-(fc2 + fc6) 











k2 


-fc3 











|fc3 


-^5 



(3.8) 



Equations (|3.6p and (|3.7I) are equivalent to the original system. The solutions 
to be constructed are obtained as fixed points of a mapping depending on pa- 
rameters. To define this mapping some manipulations of the basic equations 
p.6p and p.Tp are necessary. The equation p.6p can be rewritten as 



1. 



^ = 25fc4(a:30^(6 - 6 + ( ^^3 -kj- 



k2X2 
5X3 



(3.9) 



The first term on the right hand side can be split into a leading term and a 
remainder with the result that the whole right hand side can be written as 



25fc4(x36)=(6-e) + ^^(a;,0 



where 



F{x,0 = 25fc4(x3)*K* - a](6 - C) + 5C ( Jfc3 - fc7 - ^ 

\5 5x3 



(3.10) 



(3.11) 
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The equation (13.71) can be solved by variation of constants to give 

= exp(A'ft)S(0) + / c^p{M{t ~ s))R[x,^]{s)ds. (3.12) 



The function R[x,^]{s) depends on the functions x{s) and ^(s). 

The matrix M is of the type considered in Appendix 2. Its determinant is 
kik^k^ (fee ~ 1^2)- Thus when k2 > Sfcg it has exactly one positive eigenvalue 
a and all its other eigenvalues have real parts less than a. Moreover, there is 
an eigenvector V with eigenvalue a all of whose components are positive. Thus 

exp{Mt)x{0) = Mie"V + ^p{t; ^2, M3, M4) (3.13) 

where the components of the function ip are of the form 

4 

^i{t\ ^Ji2, ^J■3, ^J"^) = ^A^jW^^je'^^'^i^ifci''- (3.14) 

j=2 k>0 

Here the f3j are the eigenvalues of M other than a, the Wij are the components of 
the eigenvectors other than V and the Vjk are constants which are only non-zero 
for fc > if the eigenvalue f3j has multiplicity greater than one. The function ip 
depends linearly on the parameters /i2, A's and ^14 and satisfies 'ip{t; /i2, /^s, /i4) = 
Q^g(a-is)t-j as t — ?> 00 for some e > 0. 
From the formula for x we get 

a;3(t) = ^ie"*V3 + e3-?A(t;^2,M3,M4)+ / 63 • exp(Af(t - s))i?[x, ^](s)ds. (3.15) 

This can be used to split the quantity X3 in the first term of (|3.10p into a part 
containing the leading term in (|3.15p and a remainder term. The result is 

^ = 25fc4(/iil46e"*)^(6 - e) + G[x,^,t] (3.16) 

where 

G[x,tt]=25k4a [4 - (Mi^ae"*)^] (6-0 + ^0- (3-17) 

The idea is that G contains all the contributions to (13.61) which can be considered 
small. 

Equation p.l6p can also be treated by variation of constants. Let 

m = exp f-^[{f,,V,^^e"')l - l{\ . (3.18) 



Then 



m^^s+Vom+ I ^^^G[x,ts]ds (3.19) 
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with an arbitrary constant 770. Finally p.l2p can be rewritten as 



x{t) = tiie°'W + ^p{t;ti2,f^3,f^i)+ eMM{t- s))R[x,(]{s)ds. (3.20) 

Jo 

The integral equations p.l9p and (I3.20p are those which are used for the fixed 
point argument. 

Let X be the set of continuous functions {x{t),£^{t)) defined on the interval 
[0, 00) which satisfy the inequalities 



<6t^i, \m-is\e^ <K (3.21) 



for positive constants /ii, S and K which are restricted by some additional 
conditions in the fixed point argument. Denote the right hand sides of equations 
(Pl^ and ([5:^ by Ti{^,x) and 7^(^,x), respectively and let T = (Ti,7^). 
Then ^ and x solve p.l9p and (|3.20p if and only if (^, x) is a fixed point of T. 
In order to ensure that the mapping T is well-defined on X it suffices to assume 
that 5 < ^ miuj Vi and that -ft' < ^Cs since these inequalities imply the positivity 
of X and ^. The aim is to show that for a suitable choice of the constants /Xi, 
A*2, M3i VOi S and K this rule defines a mapping T from X to itself. 

For convenience let T{^,x) = (C-2/)- The quantities C and y should now be 
estimated under the assumptions p.2ip . The quantity C ~ is estimated first. 
It can be written as a sum Qi where the individual terms are defined as 

follows. Qi denotes the second term on the right hand side of (|3.19p . Q2 denotes 
the contribution to the right hand side of (|3.19p coming from the first term on 
the right hand side of p.l7p . The contribution to the right hand side of (|3.19p 
coming from the second term on the right hand side of p. 171) is Q3 + Q4 + Q5. 
The three summands in this last expression come from three summands in F. 
The expression is the contribution coming from the first term in (I3.1ip . Q4 
is the contribution from the first two terms in the bracket on the right hand 
side of (|3.1ip and Q5 is the contribution from the third term in that bracket. 
As a first step towards estimating ( ^ S,s consider the following identity which 
holds for any positive constants A and 7. 

— (7^^-16--^* exp(yle'^*)) = exp(Ae'^')(l - A-^e-^*). (3.22) 
at 

If ^ > 2 then the second factor on the right hand side of this equation is no 
smaller than one half for any t > 0. Hence integrating this relation from to i 
gives 



Jo 



(3.23) 



Thus 



\c~is\<m 
2 

25k4 



* 1 



Jo *(s) . 

< 7T^(m^36)"^l|G||Looe-^* + \vo\m- (3.24) 
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The integral has been estimated using p.23p with A — ( ^4^'' ) (Aii^^aCs)^ and 
7 = Assuming /ii sufficiently large ensures that the lower bound on A 
required for (|3.23p to hold is satisfied. Next G is estimated. 

||G|U= < 20ha K{fiiV3)^Vf^ sup{l,{V3 ~ S)-^6 + \\F\\l^. (3.25) 

Choosing S small enough allows Q2 to be bounded by -|^e~t"*. The first term 
in the expression for F can be bounded by 

Cfc4Aife^"*(F3 + <5)s(e-6)' < Ch^l!K\V3+S)ie-^' (3.26) 

for a numerical constant C. Choosing K sufficiently small allows Q3 to be 
bounded by -|^e~5"*. In the second term in F the first two contributions can 
be bounded by 2^5. (fcs + bkj). Choose fii large enough so that 

^(/^i^36)^'6(fc3 + 5kj) < K. (3.27) 

Then it follows that (54 is no greater than -|^e~i"'. The third contribution to 



F can be bounded by 2^sfc2 



Now 



To get a lower bound for the denominator in this expression it is assumed that 
6 < i^a. Then 

< \; (3.29) 



If /ii is chosen large enough then Q5 is bounded by ^^e 5"*. To control the 
quantity Qi assume that fii is so large that 

125fc4,, ,4 4q, , 4a , , 

^[(W^3e«)^e-* ^ 1] > —t. (3.30) 

for all i > 0. Then |77o|<i>(t) < |77o|e~^*. Choose r]o to be no larger in modulus 
than ^ . Then combining the estimates shows that for K and 6 sufficiently small 
and /ii sufficiently large the second defining inequality of the set X is satisfied 
by C- 

Next y — fiie°'*V is estimated. The function ip can be bounded by an ex- 
pression of the form Ce*^""*^^* where the constant C depends only on the matrix 
M and e is a positive constant. Since "0 is linear in the parameters /X2, /X3 and 
/X4 the constant C can be made as small as desired by making these parameters 
small. Thus it can be ensured that e~"*|V'(i)| < The integral term in 

(|3.20l) can be estimated by an expression of the form C/ii/sTe"'. This can be 
made as small as desired compared to <5/iie"* by choosing K small. Thus the 
first defining inequality of the set X is satisfied by y. 
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It follows that for a suitable choice of the parameters T maps X into itself. 
Let jp[t) = ^ie"*y and £,^{t) — ^s- Define a sequence recursively by 

(^"'"'"^,5;"+^) ~ a;")). The sequences a;" and ^" are uniformly bounded 

on compact subsets. It then follows from the definition of T that their time 
derivatives are uniformly bounded on compact subsets. By the Arzela-Ascoli 
theorem there is a subsequence (^"'' , x"'' ) which converges uniformly to a limit 
on compact subsets. It is possible to pass to the limit in the integral equations 
defining the iteration to see that the limit is a fixed point of T and hence the 
desired solution. The solution is uniquely determined by the parameters Hi and 
770. The mapping from these parameters to the initial data for the solution at 
t = is a diffeomorphism onto its image. Thus the set of solutions constructed 
in this way corresponds to an open set of initial data. This completes the proof 
of the theorem. 



4 Michaelis-Menten via mass action kinetics 

Next the system (10) of [9] is considered. For the convenience of the reader we 



reproduce the essential equations here: 

iRuBP = fclsa^RuSPEa " fcia^RuBpa^Ei + fc2a;RuBPEi, (4.1 

iRuBPEi = fcia;RuBpa;Ei - (^2 + fc3)a;RuBPEi, (4.2 
ipGA = SfcaXRuBPEi - fc4a;pGAa;E2 + fcsa^PGAEa 

-fci62;pGAa;E6 + fcl72;pGAE6, (4.3 

ipGAEa = fc4a;pGAa;E2 - (^5 + fc6)a;pGAE2, (4.4 

ioPGA = fcea^PGAEa - fc7a;DPGAa;E3 + fc8a;DPGAE3, (4.5 

iDPGAEa = fc7a;DPGA2;E3 " (^8 + fc9)a;DPGAE3, (4.6 

iGAP = fc9a;DPGAE3 - 5fcioa;GApa;E4 + SfcnXQAPEi 

-fci9a;GApa;E7 + fc2oa;GAPE7, (4.7 

iGAPE4 = fcloa^GAP^^Ei - (fell + fcl2)a;GAPE4, (4.8 

iRu5P = -fci32;Ru5pa;E5 + fci4a;Ru5PE5 + 3fci2a;GAPE4, (4.9 

iRuSPEs = fci3a;Ru5pa;E5 - (fcl4 + fcl5)a;Ru5PE5, (4.10 

ipGAEe = fci6a;pGAa;E6 - (fci7 + fci8)a;pGAE6, (4.11. 

^gapEt fci9a;GApa;E7 - (^20 + fc2i)a;GAPE7- (4-12 



The equations for the concentrations of the free enzymes have been omitted 
since they can easily be reconstructed. As explained below, this could be made 
into a closed system by using the conservation of the total quantity of each 
enzyme. This has not been done so as to prevent the equations becoming even 
longer. The kinetics is called Michaelis-Menten represented in terms of mass 
action in [9] and is called the system MM-MA here. The unknowns are of three 
types. There are the concentrations of free substrates, which are denoted by 
the same variables as in the system MA. There are the concentrations xe„ of 
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the seven enzymes corresponding to the seven reactions in the system. FinaUy, 
there are the concentrations of the complexes formed when the enzymes bind 
to their substrates. The complex formed when the substrate Ai binds to the 
enzyme is denoted in what follows by AiEa- Since in some reactions r{a) > 1 
molecules of the substrate bind to the enzyme this might be denoted instead 
by A^^'^^Ea. Since, however, the exponent r{a) is uniquely determined by a 
we chose the shorter notation to prevent certain formulae becoming even more 
cumbersome than they already are. It should be warned that the word 'complex' 
is used in two different ways in the literature on reaction networks. The first 
meaning is the one just introduced. The other is a formal linear combination 
of the chemical species which is on the left or right hand side of a reaction. 
To distinguish these two concepts in what follows they will be referred to as 
'enzyme-substrate complex' and 'reaction complex' respectively. 

The total concentration of each enzyme (free plus bound) is a conserved 
quantity and is denoted by pa- It follows as for the system MA that the set 
S is invariant. The question, whether solutions starting in S can have w-limit 
points on the boundary is a little more complicated than for the system MA. 
Note first that there is a five-dimensional set Ai of stationary solutions in S 
defined by setting the concentrations of all enzymes to zero together with those 
of the corresponding complexes. This is the set where all pa are zero. This 
set cannot contain any w-limit point of a solution with positive initial data, 
since for a solution of that type the pa are positive. The conservation of the pa 
defines invariant affine subspaces of S of codimension seven. It is elementary 
to show that the stoichiometric matrix has rank twelve so that these subspaces 
are the stoichiometric compatibility classes. Call one of these subspaces Se- 
Another set of stationary solutions A2 , of dimension seven, is defined by setting 
the concentration of all substrates and all enzyme-substrate complexes to zero. 
Any subspace Se intersects this set of stationary solutions in precisely one 
point. Consider a solution x{t) which is positive and which has an w-limit point 
X* on the boundary of S. The solution y{t) passing through x* lies entirely 
in the boundary of S. If the concentration of any free enzyme Ea vanishes 
at X* then, by the conservation laws, the concentration of the corresponding 
substrate-enzyme complex is non-zero. It follows that the time derivative of 
the concentration of E^ is positive, a contradiction. Thus it can be concluded 
that the concentrations of all free enzymes are non-zero at x*. Suppose that 
X* ^ A2. If all substrates had zero concentration at x* then at least one enzyme- 
substrate complex would have to be non-zero and the evolution equation for that 
substrate would give a contradiction. Thus at least one substrate must have non- 
zero concentration. Then the evolution equation for a complex of that substrate 
with any enzyme implies that the concentration of that complex must be non- 
zero. Since x* is on the boundary of S it is not possible that the concentrations 
of all substrates are non-zero. Prom this point on it is possible to argue as in 
the corresponding proof for the system MA to obtain a contradiction. It can 
be concluded that any w-limit point of a solution starting in S must either be 
a point of S or belong to A2. The conservation laws show that the w-limit set 
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contains at most one point of ^2- 

The system (|4.ip - (l4.12p has the property that all the variables representing 
the concentrations of enzymes or enzyme-substrate complexes are bounded due 
to the conservation laws for the quantities pa- On the other hand, the quantities 
on the right hand sides of the evolution equations for all other variables are 
all either non-positive or linear in the quantities which are not known to be 
bounded. It follows from this that all solutions exist globally in the future. 
Summing up; 

Proposition 2 A solution of the system (|il1) - (|4.12l) (the system MM-MA) with 
positive initial data exists globally to the future, remains positive and has no 
w-limit points on the boundary of S except possibly a single point of the set A2 . 
The conclusions listed for the system MM-MA in this proposition also hold for 
the analogous system MM-MAZ defined using the stoichiometric coefficients of 
|22j and can be proved in the same way. 

For 1 = 1,2, 3, 4, 5 let Xi be the sum of the concentration of the free substrate 
i and its concentrations within its complexes with enzymes. Note that here it 
is necessary to take into account that in general the complex contains several 
molecules of the substrate. Then A2 is the subset of S where all Xi vanish. 
These quantities satisfy the evolution equations 



dxi 
~dt 

dX2 

~dt 

dxs 

dt 

dX4 

~dt 

dX5 

~dt 



^isa^RuSPEs - fcaXRuBPEi , (4-13) 

2fc3XRuBPEi - fcea^PGAEs - fci8a;pGAE6, (4-14) 

fc6a;pGAE2 - fc9a;DPGAE3 , (4-15) 

fcga^DPGAEs - 5fci2XGAPE4 " fc2ia;GAPE7, (4.16) 

3fci2a;GAPE4 - fci52;Ru5PE5- (4-17) 



Let Li be the quantity obtained by replacing Xi by Xi in the expression for the 
function Li introduced for the system MA. Then 

dZi 1 /, 1, \ 3, „, 

= I KlS^^PGAEg - ^K6a;pGAE2 I ^ ^K212:gAPE7- (4-18) 

This shows that Li is a Lyapunov function on the region where the quantity in 
brackets in (|4.18p is non-negative. It will now be shown that this can be used 
to prove that certain solutions tend to zero as t — >■ 00. 

Proposition 3 A solution of the system MM-MA (the system (|4ll - (|4.12p ') 
with fci7 -f fci8 < fcs -|- fcg which satisfies the inequalities (|4.19p and k/^kQp2 < 
kiekisipe — 2Li(0)) converges to a point of A2 a.s t ^ 00. 

Proof For a solution satisfying the assumptions of the proposition the quantity 
in brackets in (I4.18P is initially positive, i.e. 

fc6a;pGAE2(0) < 5fci8a;pGAE6(0). (4.19) 
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The evolution equations for the concentrations occurring in this inequahty are 



da:pGAE2 



fc4a;pGAa;E2 - (^5 + fc6)a;pGAE; 



(4.20) 



dt 



rfa^PGAEs 



fci6a;pGAa;E6 - {kn + A:i8)xpgae6- 



(4.21) 



dt 



Let be supremum of times for which fc6a;pGAE2(i) < 5fci8a::pGAE6(i) holds 
on the interval [0,i,). Using the fact that kir + fcis < ^5 + ^6 the sum of 
the contributions of the second terms on the right hand sides of the evolution 
equations for a;pGAE2 and a;pGAE6 to the evolution equation for 5fci8a;pGAE6 ~ 
fc6a^PGAE2 is positive when t — t^,. Now fc4XpGAa;E2 < k4XpGAP2 and 

A^iea^PGAa^Eg = kieXpGAiPe - a;pGAE6) > fci6a;pGA(p6 - 2Li(0)). (4.22) 

Thus due to the inequality k^kQP2 < fcie^islpe ^ 2Z/i(0)) the assumption that 
tf is finite leads to a contradiction. In addition it can be seen that in this case 
any cj-limit point of the solution must satisfy xpga = and hence belong to 
A2. This gives the conclusion of the proposition. 

Next it will be shown that the system MM-MA has solutions for which 
the concentrations of the substrates tend to infinity at late times. To do this 
it is most economical to do the calculations in the framework of a class of 
reaction networks wider than those describing the Calvin cycle. Consider a 
system of chemical reactions as defined by sets of species, reaction complexes 
and reactions. This will be called the basic reaction network. It is possible 
to build a new network by replacing each reaction in the basic network by a 
Michaelis-Menten scheme containing a substrate (the species from the basic 
network), an enzyme and a substrate-enzyme complex. Applying mass action 
kinetics to the extended network gives 'Michaelis-Menten expressed in terms 
of mass action' kinetics or, for short, MM-MA kinetics. In this way starting 
from the basic network we get a system of ordinary differential equations called 
the MM-MA system. It contains reaction constants for each reaction in the 
extended system as parameters. Call the substrates Ai and the enzymes for 
some indices i and a. The complex formed when these bind to each other is 
denoted by AiEa- 

Some restrictions will now be made on the basic set of chemical reactions. 

1. Each reaction complex in the basic network contains only one species 

2. The set of substrates and the set of enzymes are disjoint. This rules out 
the MAP kinase cascade 12!. 

3. Each enzyme catalyses only one reaction. This rules out systems with 
enzyme sharing such as the multiple futile cycle |20) . 

When there are n species and r reactions in the basic network then the number of 
species in the corresponding MM-MA system \s n + 2r. There are n substrates, 
r free enzymes and r substrate-enzyme complexes. There are 3r reactions. 
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In the motivating example for this work the basic system is defined by the 
equations (|2.ip - (|2.5p describing the Calvin cycle. It has five species. There are 
seven reactions and so there are nineteen species in the corresponding MM-MA 
system. When a system is so big it is not very efficient to write it out explicitly 
when analysing it. It can be more useful to treat it as an example of a class 
of systems characterized by some particular structural properties. This is the 
motivation for considering a more general class of systems here. Note that the 
first restriction above rules out the more detailed models of the Calvin cycle 
given in [T7] and [TH] . It also rules out the homogeneous case of the model with 
diffusion considered in [S] . It will be seen in Section [5] that in fact all solutions 
of the latter system are bounded. 

The main theme of what follows is solutions of an MM-MA system in which 
the concentrations of all substrates tend to infinity as i — ?> cx). In fact they all 
tend to infinity linearly in time. In the solutions of interest here the concen- 
tration of each free enzyme tends to zero as t ^ oo and almost all the enzyme 
becomes bound to the substrate at late times. A class of networks are con- 
sidered which are called autocatalytic. It is shown that for MM-MA systems 
arising from networks satisfying this additional property, which is defined later, 
there are large classes of solutions of the type just described. They are referred 
to here as runaway solutions. 

The MM-MA system can be written as a set of evolution equations for the 
substrates, the substrate-enzyme complexes and the free enzymes. The right 
hand sides of the equations of the second and third types for a given choice of 
enzyme differ only by an overall sign. Adding them gives a conservation law 
for the total amount of enzyme pa — a^AiE„ + xe^, ■ The conservation laws can 
be used to eliminate the concentrations of the free enzymes from the evolution 
equations for the substrates and the substrate-enzyme complexes. The evolution 
equations for the free enzymes can be discarded. This leads to the system 



dt 

a:i{a)—m 

+ ^ r(a)C^ (a)xA,^^^E^ + ^ s(a)r(a)xAi,„)E„ , (4.23) 

a:i{a)—7n oc:f{a)—m 

^^^^^ - ^^+(")^A.2, (P- - ^A,„,E J - (C-(a) + r(a)):rA,„,i44-.24) 

Here C+(a), C~{a) and r(a) are the reaction constants for the reactions in- 
volving the enzyme Ea- The numbers r(a) and s{a) are the stoichiometric 
coefficients of the reaction catalysed by Ea, referred to for brevity as the reac- 
tion a. The number of molecules of substrate entering the reaction is r{a) and 
the number of molecules of product which result is denoted by s{a). In fact we 
allow r{a) and s{a) to be any real numbers satisfying the condition r(a) > 1. 
This inequality ensures that the coefficients in the system of ODE are C^. i{a) 
is the index labelling the substrate entering the reaction a and /(a) is the in- 
dex labelling the product of that reaction. The full MM-MA system consists 
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of (|4.23p . (I4.24P and evolution equations for the concentrations of the free en- 
zymes. Equations (|4.23p and (|4.24l) are equivalent to the full MM-MA system 
in the following sense. If a solution of the full MM-MA system is given then 
the conserved quantities pa can be computed. Then the concentrations of the 
substrates and the substrate-enzyme complexes satisfy (|4.23p and (|4.24p with 
those values of the pa ■ Conversely, suppose that a solution of (|4.23p and (|4.24p 
is given with certain values of the pa and that the concentration of a:Ai(„)Ec is 
always less than pa. Then defining the concentrations of the free enzymes by 
xs^ — Pa — 2^Ai(„)EQ gives a solution of the full MM-MA system. The following 
linear combination of equations (|4.23p and (|4.24p will be useful later. 




^ r(a)a;A,(„)E„ 

a:i{a.)—m 

r(a)r(a)a;Ai(„,E„ - ^ s(a)r(a)a;Ai(„,E„ ■ (4.25) 



In order to investigate when the MM-MA system admits runaway solutions 
a first step is to look for consistent leading order asymptotics. This is done 
using the following ansatz. 

XA^ - e-mt + (4.26) 
:rE„ = Vo.t-''^"'^ + .... (4.27) 

For consistency XA^^-^'Ea = Pa ~ + . . •. These relations and their for- 

mal time derivatives are now inserted into the evolution equations. Comparing 
coefficients results in the equations 

Om^- r{a)C+{a)e\[2]7ja 

OL:i{a.)—m 

+ Y r{a)C-{a)pa+ ^ r(a)s(a)pa, (4.28) 

a:i{a.)—7n a:f(a)—m 

= C+(a)0j„"))r;„ - (C-(«) + T{a))pa. (4.29) 

Substituting the second equation into the first (or comparing coefficients in 
(1123)) gives 

9^ = ~ Y Ka)r(a)Pa+ J2 5(a)r(a)pa. (4.30) 

a:i(a)—m a:f(a)—m 

Since the 9m are positive this implies a linear system of inequalities for the quan- 
tities Pa ■ If these inequalities admit non-trivial solutions then the network is said 
to be autocatalytic. For a general network it is not easy to determine whether it 
is autocatalytic. The network of [9j modelling the Calvin cycle is easily shown 
to be autocatalytic. The network obtained by replacing the stoichiometric co- 
efficients used in ^ by those used in (55] can be checked to be autocatalytic by 
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an almost identical computation. When a network is autocatalytic and the pa 
satisfy suitable inequalities then the constants 9^ are determined by equation 
(|4.30p and the constants rja are determined by equation (14.291) . 

In order to prove the existence of runaway solutions for autocatalytic MM- 
MA systems it is convenient to introduce new variables adapted to the expected 
asymptotics. Define 

XA^{t) = Z^{t){t + R), (4.31) 
a;E„ =Ca(t)(t + i?)-''("^ (4.32) 

Then the solutions to be constructed should satisfy Z„i{t) Om and C,a{t) Va 
as f cxD. The parameter R > 1 has been introduced to ensure that the leading 
terms in the quantities which tend to zero are already small for t = 0. Rewriting 
the evolution equations in terms of the new variables leads to the system 

a:i{a)—'m 

+ ^ r{a)C^{a)pa+ ^ s{a)r{a)pa - F^{Ca), 

a:i{a)—m a:f{a.)—m 

^ + C+{a){t + i?)'^(")zi";Co - r{aKt + R)-\^ 

+ (i + RY'^''\C-{a) + r(a))pa - iC-{a) + r(a))Ca (4.33) 

where 

i^m(Ca)- r{a)C~{a)Ut + R)-'-^"'> 

Q:z(a)— m 

+ J2 sia)T{a)Ut + Rr'^^^l (4.34) 

a:f{a)—m 

The main result is 

Theorem 3 Let an autocatalytic reaction network be given. Then the corre- 
sponding MM-MA system can be written in the form (|4.33l) - (|4.34l) depending 
on a parameter R. Fix the values of the reaction constants. Define parameters 
9m and rja by equations (|4.30l) and (|4.29l) and suppose that the 9„i are positive. 
Then there exist positive constants K, Rq and Sq such that if i? > Rq and 



E 



Zm{0) - 9m\ + Y IC"(0) - Vcl < So (4.35) 



then 



for alH > and 



Y \Z,n{t) - 9,n\ + Y -Va\< KSo (4.36) 



lirn ( ^ |Z,„(i) -9m\+Y ~ vM ^ 0- (4.37) 

\ m a / 
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To prove this theorem the first step is to rewrite the evolution equation for (a 
as an integral equation using variation of constants. Define 



^a{s,t) = exp 



(4.38) 



Here the fact that 4'^ depends on Zj^^j) has not been made expHcit in the 
notation. Then 

CaW = Ca(0)*a(0,i)+ / (s, t) (C" (a) + r(a) )p„ (s + i?)'-(")ds 



+ / ^c.is,t)[r{a){s + R)-\^{s)-iC-{a)+T{a))Ca.is)]ds. (4.39) 



The second term on the right hand side of this equation can be transformed 
using the identity 



1 



C+(a)(Z,(„)(s))-(") ds 
and integration by parts. The result is 
(C-(a)+r(a))p„ 



(vl/„(s,t)) = (s + i?)'-(")vl/„(s,t) (4.40) 



Ut) 



C+{a){Z,^^){t)Yi-) 

+ / *a(s,t) 



Ca(0) 



* r{a){C-{a)+T(a))po,dZ,(^){s) 



C+(a)(Z,(„)(s))-(")+i 



(C-(a) + r(a))p„ 
C+(a)(Z,(„)(0))'-(") 

ds + .. 



ds 



*a(0,t) 

• (4.41) 



where the last term in (I4.39P has not been written out again. 

Proof of Theorem 3 In this proof it is assumed that K and Rq are greater 

than one. For positive constants K and 5 define 



sup<^t >0:^|Z,„(t) 



(4.42) 



Suppose that (5o < \K^^ min{0,„}. This implies that the inequality Zm > 
holds on [0,^*]. The time derivatives of the quantities Zm can be bounded by 
a constant K depending only on the parameters in the system. To obtain esti- 
mates for the following auxiliary estimate is useful. For a positive constant 
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exp 



(r(a) + l) / {u + Ry'^"Uu 



ds 



1 d 



1 

a{r{a) + l){t + Ry(°^) 
a(r(a) + l)(s + i?)'^(") 



a(r(Q) + l)(s + i?)'-( 
Gxp[ai?''(")+i - a{t + 
a(r(a) + 

exp[a(s + - a{t + Ry^°'^+^]ds. (4.43) 



Choosing Rq large enough ensures that the first factor in the last integral is 
smaller than i. Thus the integral term can be absorbed into the left hand side 
of the inequality. Discarding a term with a good sign gives 



exp 



a{r{a ) + l) / {u + Ry'^°'Uu 



ds < 



a{r{a) + + ' 
Making a suitable choice of the constant a leads to the inequality 



[ 'i/a{s,t)ds < K{t + Ry 
Jo 



(4.44) 



(4.45) 



Putting this information into (|4.4ip gives 



{C-{a)+T{a))pa 



C+{a){Z,^^){t)yi-) 



< 



Co(o) 



(C-(a)+r(a))po 



-K{t+R)' 



C+(a)(Z,(„)(0))K") 

(4.46) 

The first term can be bounded using the smallness condition on the initial data 
and the second by using (|4.45l) and choosing Rq large. It follows that for i?o 
sufficiently large and 5q sufficiently small 



(C-(a)+r(a))p„ 



C+(a)(Z,(„)(t))K-) 



< 



K5o 



(4.47) 



The evolution equation for can be rewritten in the form 



{t + R) 



~dt 



+ (Zm ^ dm) — 



- y: r{a)C^a)Z:;:^ 

a:i{a)—m 
' Fm ( CcK ) • 



(C-(Q) + r(a))p„ 



(4.48) 
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The right hand side of (|4.48l) can be bounded by K6q /2 after possibly increasing 
K and Rq. Integrating this gives an inequality of the form 

{t + R)Y, \Zrn - e„,\ < iK5o/2){t + 1) (4.49) 

m 

and hence 



KSo 



^\Zrn~Om\<^. (4.50) 

m 

Combining (|i?5g| and (I07)) gives 

J2\Ut)-V»\<^- (4.51) 

a 

It can be concluded that t* ~ oo and the first part of the theorem is proved. 
The integrand in the definition of is bounded below by a positive constant 
and thus ^'^(O,^) ^ as t oo. Combining this with (|4.4ip shows that 

C+(a)(Zi(„)(t))n") 

ast oo. It can then be concluded from (|4.48l) that ■^{{t+R){Zm—Om)) = o(l). 
Hence (t + R){Zm — 9m) = o{t) and Z^. 0m as i — >■ oo. Together with the 
information we already have this implies that — ^> as i — ?► oo and this 
completes the proof of the theorem. 

Consider now stationary solutions of MM-MA. Equation (|4.25p implies that 
the equation obtained by setting 6'm = in (|4.30p holds in the stationary case. 
This is a linear system for the substrate-enzyme complexes. It is a system of n 
equations for r unknowns. In the case of the system (|4.ip - (|4.12p there are five 
equations for seven unknowns and it is easily seen that the solution space is of 
dimension two. The equations are those obtained by setting the time derivatives 
to zero in (|4.13p . Explicitly 

a;RuBPEi = -7— a;Ru5PE5, (4.53) 

fcl8 

a^PGAEa = -, — a^RuBPEi ^a;pGAE6 , (4.54) 

fee fee 
a;DPGAE3 = 7^a;pGAE2, (4.55) 

fcg 

fcg ^21 

a;GAPE4 = 71 a^DPGAE, - —, XGAPE7, (4.56) 

5fel2 ' 5fei2 

a;Ru5PE5 = -r — a;GAPE4- (4-57) 

Suppose now that we prescribe the values of xpgaEb and xgapEt ■ It is possible 
to derive the equation 

a;GAPE4 = 7^(^i8a:^PGAE6 + fc2ia;GAPE7)- (4.58) 
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Substituting back into equations (|4.53p - (|4.55p and (|4.57p gives: 



a^RuBPEi 


3 

— 1— (fcisa^PGAEe + feia;GAPE7), 


(4.59) 


a^PGAEa 


= — (5/ci8a;pGAE6 + 6A;2ia::GAPE7), 
fee 


(4.60) 


2;dpgae3 


= — (S/cisXpGAEe + 6A;212:gAPE7), 
fcg 


(4.61) 


a;Ru5PE5 


3 

= T— (fclS^^PGAEs + feia;GAPE7)- 


(4.62) 



These can then be used to obtain expressions for the concentrations of the 
free enzymes. For the total amount of any one of the enzymes is equal to the 
amount of the free enzyme plus the amount of it bound to its substrate. Now 
these expressions will be used to extract information from the time evolution 
equations for the free enzymes. For brevity let X = xpcAEe ^^id Y = xgapEj- 
Then 



a^RuBP 



XPGA 



a^DPGA 



XGAP 



a^RuSP 



3(fc2 + fc3)(fci8X + fc2iy) 

ki{hPi - 3fci8^ - 3/c2i5^) ' 

(fc5 + fc6)(5fci8X + 6fc2iy) 

ki{kQP2 - SfcigX - &k2iY) ' 

(fc8 + fc9)(5fcl8^ + 6fc2iy) 

kjikgps - 5kisX - 6/c2iF) ' 
{kii+ki2)iki8X + k2iY)' 



kw{ki2Pi - fcis-'^ - k2iY) 

3(fci4 + fcl5)(fcl8^ + fc2iy) 
fcl3(fcl5P5 - 3fcl8-'^ - 3fc2iy) ' 



(4.63) 
(4.64) 
(4.65) 

(4.66) 
(4.67) 



The expressions obtained up to now suffice to determine all unknowns in terms 
of X, Y and the conserved quantities pa- There are, however, two further 
equations which lead to consistency conditions. These are: 



a;pGA 



a;GAP 



(fci7 + fcis)^ 
kw{p& - X) ' 

(fc20 + fc2l)y 

kig{p7-Y) ■ 



(4.68) 
(4.69) 



Note that equations (|4.63p - (|4.69p only have positive solutions under the restric- 
tions that X and Y satisfy the inequalities which ensure the positivity of the 
denominators of the right hand sides. Rearranging the equations (|4.68p and 
(|4.69p gives 



X 



Y 



kiePeXpGA 



f kis + fcigXpGA ' 

kmPrXGAP 



Cl7 



^20 + k2i + fcigXGAP 



(4.70) 

(4.71) 
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Hence 



bkisX + 6k2iY = + , _^ , ^ , (4.72) 

ki7 + fcis + Kiea^PGA K20 + fc2i + Ki9a;GAP 

and 

fcigX + k2xY = + , ^ , . , ■ (4.73) 

fci7 + + fci6a;pGA fc20 + K21 + fci9a;GAP 

Substituting these relations into (|4.64l) and (|4.66|) gives equations of the form: 

a;pGA - 51 (a^PGA , a^GAP ) = 0, (4.74) 
s^GAP - 52(a;pGA, a;GAp) = (4.75) 

for some rational functions g\ and 52. More explicitly 

, . aia;pGA + a2a;GAP + aaXpGAa^GAP , , nn\ 

5i(a::pGA,a;GAp) = 7— -T TT -T > (4.76) 

bi + 62a;pGA + Osxgap + 04a;pGAa;GAP 

/ s CiXpGA + C2a:GAP + CsXpGAa^GAP I , nn\ 

32(XY>GA,XGKY>) = , , , —, —, (4.77) 

di + d2a;pGA + daa^GAP + a4a;pGAa;GAP 

for suitable constant coefficients depending on the ki and the p^. 
Lemma 1 Any positive stationary solution of the system (j4.ip - (|4.12p defines a 
positive solution of the system (|4.74l) - (|4.75p . Conversely each positive solution 
of (|4.74p - (|4.75p with given values of pa for which the quantities X and Y defined 
by (I4.70p and (|4.7ip make the denominators in (|4.63p . (|4.65p and (|4.67p positive 
defines a positive stationary solution of the system (|4.ip - (j4.12p . 
Proof The first statement is a direct consequence of the calculations which 
have just been done. To prove the converse let (xgap, sjpga) be a solution of 
(ii7il) - (|i75)) and let X and Y be defined by and (HTT|) . Then and 

((4:69)) are satisfied. It follows from (|474)) - (|475)) that (|4J4)) and (14:66)) hold. 
Next define xrubPi a;DPGA and xruSP by ()4.63p . ()4.65l) and (14.671) respectively. 
Define the quantities by the conservation laws and the quantities xruBPEi j 
a;pGAE2, ajDPGAEa: ^^GAPEi and a^RuSPEi by (|4.58p and ()4.59l) - (|4.62l) . Now ah 
the variables in the system (|4.ip - ()4.12p have been defined and it remains to 
show that they define a stationary solution. Equations (I4.58P and ()4.59p - ()4.62p 
imply ()4.53l) - ()4.57l) . At this point it is useful to think of the system ((4A)) - ()4.12p 
as a special case of the MM-MA system introduced for a more general class of 
networks above. In that framework what has been obtained up to this point in 
the proof is a stationary solution of the equations (|4.23p . (I4.24p and (|4.25p . It 
follows from the discussion above that this set of equations is equivalent to the 
full MM-MA system and this completes the proof of the lemma. 

In elementary flux modes of this system are investigated. This is a concept 
coming from chemical reaction network theory which can sometimes be used to 
investigate the number of stationary solutions of a dynamical system coming 
from a network of chemical reactions. To describe this in more detail it is 
necessary to introduce the notion of the deficiency of a reaction network. First 



23 



note that a directed graph can be associated to any network where there is a 
vertex corresponding to each reaction complex and an arrow corresponding to 
each reaction. The connected components of this graph are called linkage classes. 
Let n be the number of reaction complexes, s the rank of the stoichiometric 
matrix and I the number of linkage classes. Then the deficiency of the network 
\s 5 = n — s~l. If the deficiency of the network (which is always non- negative) is 
equal to one and some other technical conditions are satisfied information about 
the number of stationary solutions can be obtained using what is called the 
deficiency one algorithm (DIA) |S]. There is a computer implementation of this 
algorithm which can be applied to cases where the network is not too large [5] . 
An elementary fiux mode defines a subnetwork which is always of deficiency one 
[5]. Under suitable technical conditions stationary solutions of the subnetwork 
lead to corresponding stationary solutions of the full network and this can be 
proved using the implicit function theorem. In [3] this procedure is cited to 
conclude the existence of two distinct stationary solutions of the system MM- 
MA in a given stoichiometric class. In what follows we will not say much more 
about this approach but the results obtained by using it were the starting point 
of the more direct proof of the existence of two stationary solutions given here. 

In the MM-MA model for the Calvin cycle there are two stoichiometric 
generators and each defines a subnetwork. The subnetwork can be obtained 
by setting some of the reaction constants to zero. Here we procede directly 
using certain limits for the reaction constants corresponding to the two modes. 
Consider the system obtained from the system for stationary solutions of the 
system ()4.ip - (l4.12p by setting /cig, fci7 and /cis to zero. Call it LSI. If the 
limiting values of the parameters are approached in such a way that kn/kiQ 
and kis/ki6 tend to non-zero limits qir and qis then the system varies in a way 
which is smooth up to the boundary. A similar system LS2 can be obtained by 
letting kiQ, k2o and fc2i tend to zero while fc2o/^i9 and k2i/kig tend to non-zero 
limits q2o and 521- 

Lemma 2 Consider the system LSI and suppose that kQP2 — 6/1:21^7 > 0. If 
ki2P4 — k2iP7 > then there is a unique positive solution. If ki2Pi — ^21^7 < 
then the number of solutions is zero, one or two according to whether 



4^12^4(^20 + fc2i) 



5fci9(fci2/04 - ^21/37) 



kioPi{k2a + k2i) - (fell + ki2)kigk2iP7 (4.78) 



is negative, zero or positive, respectively. 

Proof In this case the functions gi and 172 in (|4.74l) and (|4.75p only depend on 
xgap- This means that (|4.75|) is an equation for xqap alone and for a suitable 
solution of this equation a corresponding value of xpcA can be calculated. The 
explicit form of (I4.75P is 

^4 ^ (fell + ki2)kiQk2iPr 

'^^^ fclofcl2P4(fc20 + fel) + fclofcl9(fcl2P4 " fc2lP7)a;GAP ' 

When fci2P4 — k2iP7 > the right hand side of this equation is non-increasing 
and the equation (|4.79p has unique positive solution. The explicit form of (|4.74p 
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IS 

^^^^ _ 6(fc5 + k(i)kiQk2iP7XGAP ,^ 

^4^6^2(^20 + fc2i) + k4ki9{kep2 - 6/e2iP7)a;GAP ' 

If the right hand side of this is positive it defines an acceptable solution for xqap 
and the first statement of the lemma follows. Consider the case ki2Pi — ^21^7 < 
0. Rewrite the equation schematically as 

= (4.8f) 



where a > 0, /3 > and 7 < 0. If p{x) — jx^ + /3a;'* then the equation to be 
solved is p{x) — a. The function p has a unique maximum at a;, — and 

p{x^) = i/3 f l^j . Comparing this quantity with a gives the second result of 
the lemma. 

Lemma 3 Consider the system LS2 and suppose that fci2P4 — kisPe ^0. It 
has no positive solution if the quantities k^kelkn + kig,)p2 — 5(fc5 + kQ)kiekigpQ 
and —kQp2 + Sfcig/Oe are non-zero with opposite signs and exactly one positive 
solution when they are non-zero and have the same sign. The solution is given 

by 

k4,ke{kir + kis)p2 - ^{k^ + kG)kiekispe ,. 
k4:ki6(-kep2 + okispe) 

^5 _ (fell -I- fci2)fci6fci8P6a;pGA 

kwki2Pi{kn + kis) + kiQkie(ki2P4 - fci8P6)a^PGA 
The only other case where there exist positive solutions is when 

kiikn + kis) ^ {k^ + k6)kie (4.84) 

and 

kep2=5kisP6. (4.85) 

In that case a;pGA is arbitrary and there is a continuum of solutions. 
Proof This is a direct calculation. 

The parameters which are contained in 171 and (72 sue 

^4, fcs, ke, fcio, fcll, fci2, fci6, fci7, fci8, fci9, ^20, ^21, P2, /34, P6, Pi- (4.86) 

These equations do not depend on pi, p^ or p^. Thus the conditions required 
to ensure the positivity of the solutions of (I4.63p . (|4.65l) and (j4.67p can be 
guaranteed by choosing pi, p^ and p^ large enough while keeping the other 
parameters fixed. In the two limiting cases considered above these functions 
simplify. In each limiting case the functions gi and 32 depend on only one of 
the variables xpga and a;GAP ■ Thus in that case one of the equations to be solved 
involves only one of the unknowns. If it can be solved then it can be substituted 
into the other equation to get the other variable. The derivative of the mapping 
sending the unknowns to the right hand side of the two equations is invertible 
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in the two limiting cases. It follows by the implicit function theorem that for 
parameter sets close to those of the limiting cases the number of stationary 
solutions is independent of the values of the parameters. It follows in particular 
that there an open set in the space of parameters and conserved quantities for 
which this construction proves the existence of two distinct positive stationary 
solutions for given values of the parameters and conserved quantities. 

It is possible to define a system MM-MAZ with mass action via Michaelis- 
Menten kinetics starting from the data in [52] . Like the other systems considered 
up to now it has the property that S is positively invariant. In the case of 
the system MM-MAZ the quantities Xi satisfy evolution equations which are 
the same as those satisfied in the case of the system MM-MA except that the 
coefhcients 5 and 3 are replaced by 1 and |. Next stationary solutions of MM- 
MAZ will be considered. The concentrations of the substrate-enzyme complexes 
satisfy a system of linear equations similar to those in the case MM-MA. The 
relation between these linear systems can be expressed succinctly by saying 
that ki2 is replaced by ^ki2- These linear equations can be solved for the 
concentrations of the first five complexes in terms of the other two. The result is 
similar to that for MM-MA with slightly different coefficients. In the equation 
for XGAPE4 there is an extra factor of five while the equations for the other 
complexes are as before since they are independent of fci2. It is also possible 
to derive expressions for the concentrations of all substrates. These are all 
identical to the corresponding equations in the MM-MA case except for the 
equation for the concentration of GAP. This last equation is changed in two 
ways. The first is that ki2 is replaced by |^i2- The second is that the exponent 
i is replaced by one. Equations similar to (j4.74p and (|4.75l) can be derived, 
with the important difference that the fifth power in (|4.75p is replaced by the 
first power. Thus (|4.79p is replaced by a linear equation. This linear equation 
has a unique positive solution provided a certain sign condition is satisfied and 
no positive solution otherwise. Thus for the system MM-MAZ, in contrast to 
the system MM-MA, this construction does not lead to a proof of the existence 
of more than one stationary solution for any value of the parameters. Setting 
the right hand sides of the equivalents of equations (|4.64p and (|4.66p for the 
system MM-MA equal to the corresponding quantities coming from the right 
hand sides of equations (|4.68|) and (|4.69l) respectively shows that the set of 
X and Y defining stationary solution is the intersection of the zero set of two 
quadratic polynomials in X and Y. Hence by Bezout's theorem [TU] unless there 
is a continuum of solutions there are at most four. 

5 Michaelis-Menten kinetics 

Starting from the MM-MA system it is possible to obtain a simplified system, 
the Michaelis-Menten system (MM system) by passing to a singular limit. This 
will now be carried out formally. It will be convenient to describe this in the 
context of the more general system introduced in the last section. (The basic 
scheme is explained in a simpler case in Appendix A.) Let e be a positive real 
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number and define r — et, iA^^jEa — e ^^M{c,)'e^c ^'^d xe^ — e '^''Ce„- Defining 
pa = e~^Pa allows the transformation of the system to be carried out directly on 
the system (|4.23p - (|4.24p . Equation (I4.23P is identical to what it was originally 
except that the original quantities are replaced by the transformed quantities. 
The factors of e cancel. On the other hand the equation (|4.24p picks up an extra 
factor of e on the left hand side. Formally taking the limit e ^ results in the 
equation obtained by setting the time derivative in (|4.24p to zero. Solving this 
equation for iAi(<:,)E„ gives 



C+{a)xJ^\pa 

5au„iE„ = • (5.1) 

" C+(a)x:^;^\ +C-(a)+r(a) 

Substituting this back into equation (|4.23l) gives the MM system 
dxA.„ _ C'+(")Par-(a)r(a)xAj"], 



C+ ia)pas{a)T{a)x''j^^'^ 



Consider the ansatz XAn, = OmT+ ■■■ for runaway solutions of ()5.2p . Substituting 
into the equation and comparing coefficients gives the equation obtained from 
(|4.30p by replacing 6m and pa by Om and pa. Define a new variable Zm by 

a;A.„(r) =Z„(T)(r + i?). (5.3) 

Then equation (|5.2p becomes 



Par{a)T(a) 



dr 



«.M=™ 1 + iCna))-^Z-:\'^>iC-ia) + r(a))(r + i?)-K«) 

Pas{a)T{a) 

"«:/M=™ 1 + {C+{a))-^Z-:\''\c-{a) + r(a))(r + i?)-(«) ' 



(5.4) 



Theorem 4 Let an autocatalytic reaction network be given. Then the corre- 
sponding MM system can be written in the form (|5.2p depending on a parameter 
R. Fix the values of the reaction constants. The positive parameters 9m are 
determined by the reaction constants. There exist positive constants K, Rq and 
6o such that ii R> Ro and 

Y,\Zm{0)~em\<So (5.5) 
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then 

^|Z„(r)-0„,| <X,5o (5.6) 

m 

for all T > and 

lim V|Z„(t)-0,„| =0. (5.7) 

T—^OO '—^ 

m 

Proof The proof is similar to that of Theorem 3 but simpler. Define 

r* = sup J r > : ^ \Z„,[t) - K^\ < ^KSo \ ■ (5.8) 



Choose i5o small enough so that |Z„j(t)| > 9/2 on [0,t*]. It follows from (15. 4p 
that 



[r + R}^ + (Zr^ ^ e^) 



<K{t + R)-\ (5.9) 



Integrating this in time and choosing Rq large enough gives an inequality of the 
form 

(r + R){Z„, - 0;„)(t) < KSoil + tY (5.10) 

for any e > 0. This allows the bootstrap assumption to be improved if r* is finite 
and it follows that in fact t* — oo. Using (|5.9p again gives the final statement 
of the theorem. 

Stationary solutions of the MM-MA system give rise to equations similar to 
those defined by the runaway solutions. In that case the analogue of (|5.ip with- 
out tildes holds. Substituting this into the evolution equation for the substrates 
gives the equation (|5.4p without tildes. This means that any stationary solution 
of the MM-MA system defines a stationary solution of the MM system. Con- 
versely, any stationary solution of the MM system defines a stationary solution 
of the system (|4.23l) - (|4.24l) . It was already shown that any solution of the latter 
system defined a solution of the system MM-MA and if the solution of (|4.23p - 
(|4.24l) is stationary the solution of the system MM-MA is so too. Thus there 
is a one to one correspondence between stationary solutions of the MM system 
and stationary solutions of the MM-MA system for fixed values of the conserved 
quantities pa ■ To get the standard form of the Michaelis-Menten system as used 
in [22] the numerators and denominators in all the summands on the right hand 
side should be divided by (7+ (a). Hence if the MM system is given in isolation 
the reaction constants of the MM-MA system it comes from are not determined 
uniquely. Only the expressions T{a) and c+(a^^"'* determined. 

Equations for a model with Michaelis-Menten kinetics are written in Ap- 
pendix A of |22j . In fact, as has been remarked in [4], the expression for 
in [22] is not correct since it includes a dependence of the reaction on ATP, 
which does not agree with the reaction it is supposed to model, the sixth re- 
action in Table B2 of [22]. Here we consider the correct model obtained from 
the reaction network given in [22] by using Michaelis-Menten kinetics and call 
it the model MMZ in what follows. The discrepancy just mentioned only affects 
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the biological interpretation of the constants in the system studied here and 
not the general mathematical form of the equations. The variables used in [22] 
include the concentration of ATP but this is assumed to be constant and so no 
evolution equation is required for it. A system with Michaelis-Menten kinetics 
is also considered (but not written explicitly) in [5] and is called MM in what 
follows. The only difference between these two systems is that the expression 
Vi in [35] is replaced by an expression of the form ^""^^ for some positive 
constants A and B. For both systems the positive orthant S is invariant. In 
these Michaelis-Menten systems the right hand sides of the evolution equations 
are bounded functions of their arguments and so solutions exist globally to the 
future. Making use of these facts it can be shown as in the case of the systems 
MA and MAZ that there are no cj-limit points on the boundary except possibly 
the origin. To see this it is enough to use the structure of the first equation in 
Appendix A of [22]. 
In both cases 

^ = -2(^^5--^^2)-5.6 (5.11) 
where the Vj denote reaction rates. With the corrected value of we get 



1 _ Vsmaxa^PGA 1 
5 (XPGA + Kmb) 5 



V2maxa;pGA2^ATP 



(a;pGA + -fi^m2l)(a;ATP + Km22) 



(5.12) 



The positivity of this is equivalent to the inequality 



(a;pGA + -K'm5)V2maxa:ATP 

< 5(a;pGA + J'^m2l)(a;ATP + ^m22)V5„iax- (5.13) 

This can be rewritten as 

5i^r„j2i(a;ATP + ^m22)V5 max 

— i^m5V2maxa^ATP /r t a\ 

^2maxa;ATP " 5(XATP + Km22)V5 

max 

provided the denominator is non-zero. Call the right hand side of this equation 
K and suppose that K > 0. Then if a solution initially satisfies the inequality 
|Li < K it tends to zero as t ^ oo. For the values of the parameters given in 
Appendix B of [22] if is about 0.14. 



6 A model including a dynamical description of 
ATP 

In [9] a model with diffusion is considered which is given by the equations (13) 
of that reference. They define a system of reaction diffusion equations. This 
system is denoted by MAd. Setting the diffusion coefficient equal to zero (or 
restricting to spatially homogeneous solutions) gives a system of ODE different 
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from (|2.ip - (l2.5p due to the inclusion of the concentration of ATP as a dynamical 
variable. Call this ODE system MAdh. The explicit form of this system is 



dXRu 



BP 



dt 

dxpGA 

dt 

dXBPGA 

dt 

dxGAP 

dt 



= fc5a;Ru5pa:ATP - hx^uBP, (6.1) 

= 2A;iXRiiBP - fc2a:pGAa;ATP - fcea^PGA, (6.2) 

= A:2XpGAa:ATP - ^3a;DPGA, (6.3) 

k^xuPGA - 5k4x%j^p ~ krXGAP, (6.4) 



^^^^.^^ = -fc5a;Ru5P2:ATP + 3kiX%j^p, (6.5) 

—^IL = -A:2a:^pGAa:ATP - fcsa^RuSP^^ATP + kgic- xatp) (6.6) 
dt 

for a constant c. Adding equations (16.11) and (16.61) gives 

^(a;RuBP + a^ATp) = -fcia^RuBP - A:2a;pGAa;ATP + ks{c - xatp)- (6.7) 
dt 

Let m = minjfci, ks}. Then 

^(a;RuBP + a;ATp) < ~m{xKuBP + xatp) + ksc. (6.8) 

It follows that XRuBP + a; ATP can be bounded by the maximum of its initial 
value and the quantity Call this xruBP- The evolution equation for xpqa 
implies that this quantity is bounded by the maximum of its initial value and 
ipGA — 2^i^a!isz.. Similarly the quantities Xupga and xqap are bounded by 
the maximum of their initial values and 

fc2ipGAiRuBP . fcs^DPGA „\ 

a;DPGA = , xgap = , (6.9) 

Ks kr 

respectively. 

Obtaining a bound for xru5p is more complicated. Integrating the evolution 
equation for xdpga on the interval [s, t] and using the fact that a;pGA and xdpga 
are bounded leads to an inequality of the form 



^ xuPGAiOd^ - ^ (/ ^ATP(^)'^^ + ^ 



(6.10) 



for a positive constant C. Similarly it follows from the evolution equation for 
a; gap that 



^ XGAp(Orff ^ ^ (/ ^DPGA(Orfe + l) • (6.11) 

Combining these two inequalities gives 

j\GApiOd^<c(^j'^ XATp{m + ^^ ■ (6.12) 
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By variation of constants the evolution equation for xjiu5P implies 



a;Ru5p(i) = a:Ru5p(0)exp ^-fcs J XATp{s)ds 
+3^4 J exp^-fc5^ XATpiOd^^ (a;GAp(s))^ds 



(6.13) 



for a positive constant C". The first term is bounded and the second can be 
bounded by an expression of the form 



C J exp^-C"^ XGApiOd^"^ XGAp{s)ds. 



(6.14) 



The integrand in the last expression can be written in terms of the derivative 
of exp (^—C xqapIO^C^ ■ Thus this expression can be bounded by 



C 



l-exp(^-C'^ XGAp(e)de 



(6.15) 



It follows that a;R,u5P is bounded. 

Since solutions of the system MAdh are bounded they exist for all future 
times. If a solution of MAdh has an w-limit point on the boundary of 5* then 
a;ATP is positive there. It follows, using the same argument as was applied to 
the system MA, that if a solution of MAdh has an w-limit point on the bound- 
ary then all concentrations except that of ATP are zero there. The evolution 
equation for xatp then shows that xatp = c at that point. Linearizing about 
the point (0, 0, 0, 0, 0, c) shows that it is a hyperbolic sink. 

Consider the stationary solutions of MAdh, i.e. the homogeneous stationary 
solutions of the MAd model. Combining the first and fifth equations gives 
^GAP ~ (fci/3fc4)a;RuBP, just as in the system MA. In fact, as shown in for a 
stationary solution all other concentrations can be expressed in terms of xruBP • 
Note first that 

^.^^p — c - ^^l^R^BP _ ^ f klfRMBp \ ^ j-g -j^g-j 

3^8 fcg \ 3fc4 / 

For an admissible solution the right hand of this equation must be positive. If 
this is satisfied the concentrations other than xruBP can be computed. The 

result IS tnat xpga - fe^^JATP+fce ' ^dpga - ks + ka ^ x^iudP - k^xATP " 
Substituting all these relations into the evolution equation for PGA gives an 
equation for xrubp alone. It is of degree ten and hence it is not easy to extract 
information from it. 

An alternative approach is the following. Another equation which can be 
derived for stationary solutions of MAdh is 



a;GAP 



k7{k2XATP + fce 



kA{k2XATP - 5^6) 



(6.17) 



31 



On the other hand, the equation (|6.16p can be rewritten as 



a;ATP = c - -t—xqp^p - — XGAP- (6.18) 

Thus we have a set of two equations for the two quantities s^gap and xatp 
and solving these is equivalent to determining all stationary solutions of MAdh. 
Write the equation (|6.17p schematically as xqap — Fi^xatp)- The function 
Fi is decreasing on the region xatp > 5^6/^2 where it is real. Note that the 
right hand side of (|6.18l) is a decreasing function of xqap and thus this equation 
can be inverted to write it schematically as a;GAP = ^2(2:atp) for a decreasing 
function F2 defined on the interval [0, c] with i^2(c) = 0. Stationary solutions of 
MAdh are in one to one correspondence with intersections of the graphs of -Fi 
and F2. The function Fi is strictly convex since F" > 0. On the other hand, 
the function F2 is strictly concave. It follows that the two graphs can intersect 
in at most two points for any given values of the parameters. If c < 5kQ/k2 
they do not intersect at all. For fixed values of the reaction constants if c is 
sufficiently large the graphs intersect in two points. 
For the system MAdh 

= ^^6 - ■^fc2a;ATP^ XpGA ~ ^kjXGAP- (6.19) 

Note that in this system xatp is bounded by c. Hence to make Li a Lyapunov 
function it suffices to require the inequality ck2 < S/cg • Thus when this inequality 
is satisfied all solutions converge to the origin as i — )■ 00. For the system MAdh 
the function L2 with a = satisfies the same equation as it does for the system 
MA. Hence the same conclusion can be drawn about solutions which converge 
to zero. 

One remark will be made on the behaviour of solutions of systems including 
a diffusion term. In fact we can do this in any space dimension, adding diffu- 
sion terms in any subset of the equations with any choice of positive diffusion 
constants. Suppose that the domain of the spatial variables is a bounded region 
fl with smooth boundary and assume that Neumann boundary conditions hold. 
Let £1 = J^Li. Then 



( h - -fexATP ) a^PGA - ^fcra^GAP 



2 V " 5 ^^^^ '^'^ 5 



(6.20) 



When the inequality fcg — ^fc2XATP ^ holds Ci is a Lyapunov function. In 
particular, there are no stationary solutions, homogeneous or inhomogeneous, 
when cfc2 < Sfcg- 

7 Conclusions 

An important motivation for the work of [22] and [9] on models of the Calvin 
cycle was to see if photosynthesis can work in different steady states. This led to 
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the question of whether the relevant mathematical models admit more than one 
stable positive stationary solution. For related work see also [21], [4], [13] and 
[H] . It was already shown in [S] that a simple model using mass action kinetics 
(the model called MA) never admits even one solution of this type. Depending 
on the values of the reaction constants, either there is no positive stationary 
solution at all or if there is it is unstable. This suggests that this is not a 
very good model. Trying to understand more about what actually happens in 
this model leads to the mathematical question of what the solutions departing 
from a small neighbourhood of the stationary solution actually do or, more 
generally, what the long-time behaviour of solutions is. Theorem 1 of this paper 
provides partial answers to this question. When there is no stationary solution 
the concentrations of all substrates tend to zero. When there is a stationary 
solution some solutions have the property that all concentrations tend to zero 
while others have the property that all concentrations tend to infinity (runaway 
solutions) . The latter alternative seems to be a further undesirable property of 
the model. 

The model MM-MA, which has a much larger number of unknowns, does 
admit more than one stationary solution for certain values of the reaction con- 
stants, as was shown in [5], and numerical results indicate that one of them 
is stable. In this paper this existence result was reproduced by a more di- 
rect method. The model MM has the same stationary solutions as the model 
MM-MA and thus also admits two stationary solutions for certain choices of 
parameters. Applying the same method to the related system MMZ with the 
stoichiometric coefficients taken from [22] does not reveal the presence of mul- 
tiple stationary solutions and it may be that in that case there are none. From 
this point of view the models MM-MA and MM look better that the model 
MA but in fact, as shown in Theorems 3 and 4 of this paper, both the models 
MM-MA and MM exhibit runaway solutions. An intuitive explanation for the 
existence of these solutions is that in all these models the concentration of ATP, 
which is the energy source for the reactions, is not modelled explicitly. Instead 
ATP is implicitly assumed to be plentiful and thus present at a constant level. 
In [3] another model is considered where diffusion is taken into account and the 
concentration of ATP is modelled dynamically. Restricting to spatially homoge- 
neous solutions leads to a system of ODE called MAdh. Interestingly, we were 
able to show here that all solutions of MAdh remain bounded, so that there 
are no runaway solutions in that model. Although heuristically plausible this is 
subtle to prove. For suitable values of the parameters this system also admits 
two positive stationary solutions. 

There seem to be few general results available on the boundedness of so- 
lutions of systems of ODE arising from chemical reaction networks with mass 
action kinetics. One theorem says, using the language of Chemical Reaction 
Network Theory, that in a mass action system which is weakly reversible and 
has only one linkage class all solutions are bounded [2 . Neither of the main hy- 
potheses of that result hold for any of the systems considered in this paper but 
perhaps some of the techniques used there might be adapted to give information 
about models for the Calvin cycle. For weakly reversible systems it might be 
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possible to prove that solutions do not have cj-limit points on the boundary 
of the positive orthant. Information about this and relevant references can be 
found in [3]. 

It is of interest to compare the conditions which allow the conclusions of the 
theorems in this paper to be obtained with values of the parameters which are 
biologically reasonable. To do this we start from the biological data collected in 
Appendix B of (22] . In that paper values are given relating to Michaelis-Menten 
kinetics. Assuming that Michaelis-Menten kinetics goes over into mass action 
kinetics when the concentration of the substrate is small compared to that of 
the enzyme it is possible to get values for the reaction constant in equations 
dlUl-lM]). The results are 

max max max max 3max^ATP \ 

\ Kml ' Km2lKm22^ Km3 ' Km4 ' Km5 ' ^m6 ' Km,13lKjnl32 J 

= (3.78, 125, 10.1, 9.63, 4, 0.02, 4.52). (7.1) 

With these values of the reaction constants the ratio ^ which plays a key role 
in Theorem 1 is equal to 1250. Thus these values are well within the regime 
where a positive stationary solution exists. 

Is it true that for the system MA with k2 > Sfcg every solution cither tends 
to infinity, the origin or the positive stationary solution? This is not known but 
since (|2.1l) - (|2.5p is what is called a monotone cyclic feedback system it follows 
from the main theorem of |15| that bounded solutions have w-limit sets which 
are no worse that those of a two-dimensional system. Using similar ideas it can 
be shown that almost all bounded solutions converge to the stationary solution. 
The system (|2.ip - (|2.5l) satisfies > for i j and is thus cooperative. It 
is also irreducible in the sense that no non-trivial coordinate hyperplane is left 
invariant by the Jacobian at any point. Using this the result on convergence of 
all bounded solutions except for those whose initial conditions belong to a set 
of measure zero follows from a theorem of Hirsch [TT] . 

In this paper a variety of different results have been proved about the dy- 
namics of solutions of mathematical models for the Calvin cycle. A number 
of interesting open questions remain to be investigated. It would be desirable 
to have rigorous results on stability of the stationary solutions and the (non)- 
existence of periodic solutions. Inhomogeneous solutions of the system with 
diffusion should be investigated following the initial work in 9 . It would be 
good to extend the results of this paper to more general models of the Calvin cy- 
cle such as those of |T7] and jT^. Finally, the basic motivating question remains: 
are there mathematical models of the Calvin cycle where it can be proved that 
there are at least two homogeneous stable positive stationary solutions? 

A Michaelis-Menten theory 

Consider a simple reaction which converts one molecule of the species S (the 
substrate) to one molecule of the substance P (the product). With mass action 
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kinetics this leads to the equations xs = —kxs and xp = kxs- Suppose now 
that this reaction is catalysed by an enzyme E. Then there is a reaction with 
reaction constant fci in which the substrate combines with the enzyme to form 
a complex SE. The reaction constant for the process of dissociation of SE into 
S and E will be denoted by Finally there is the reaction in which the 

complex gives rise to the product with reaction constant k2 while setting free 
the enzyme. This gives rise to the system 

Xs = -kixsXE + k-ixsE, (A.l) 

xsE = kiXsXE - {k-1 + k2)xsE, (A. 2) 

Xe = -kixsXE + (fc-i + k2)xsE, (A.3) 

Xp = k2XsE- (A.4) 

The first three of these equations form a closed system and thus it is natural 
to analyse it first and use the last equation to determine the evolution of the 
concentration of the product afterwards, if desired. The above system is the 
MM- MA version of the original simple reaction. The Michaelis-Menten kinetics 
will now be derived on a heuristic level. Note first that the quantity xse + xe 
is conserved. Call it Eq. Substituting the; relation xe = Eq — xse into the first 
two evolution equations gives a closed system for xs and xse' 

xs = -kiEoxs + (fc_i + kixs)xsE, (A.5) 
XSE = kiEoxs - {k-i + kixs + k2)xsE- (A.6) 

Now introduce r = et, xse = e~^xsE and Eq = e~^Eo for a constant e. This 
gives 

x's = -kiEoXs + {k-i + kixs)xsE, (A.7) 
ex'sE = kiEoxs - {k-i + kixs + k2)xsE (A.8) 

where the primes denote derivative with respect to r. In the last system it 
is possibly to formally pass to the limit e — J- 0, corresponding to a very small 
amount of enzyme. In the limit the second equation reduces to the algebraic 
equation 

kiEoXs „^ 

XSE = T — T —j-- (A-9) 

k-i + kixs + k2 

Substituting this back into the evolution equation for xs and gives the effective 
Michaelis-Menten equation 



kixs + k-i + k2 



It can then be computed that in this set-up x'p = —x'g. 

This type of discussion is quite standard and the reason it is reproduced here 
is to illuminate the relations between the three types of kinetics (MA, MM-MA 
and MM) by an explanation of the simplest example. There is a one to one 
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correspondence between stationary solutions of the systems MM-MA and MM, 
as will now be shown. If a stationary solution {xs, xse) of the system MM-MA 
is given then the equation (jA.9p is satisfied. Hence the equation (|A.10p holds 
and a stationary solution of the system MM is obtained. Conversely, suppose 
that a solution {xs,xse) is given. Then a stationary solution of the system 
(jA.Sp is obtained. Defining xe = Eq — xse for a, positive constant Eq completes 
it to a stationary solution of the system MM-MA. 

B A special class of matrices 

This appendix is concerned with the algebraic properties of some matrices of a 
special form which appear in this paper. Let A he an n x n matrix with entries 
aij. Suppose that an < for each i, that > for j = i — 1 mod n and that 
Uij = otherwise. Suppose further that ( — 1)"+^ det ^ > 0. The matrix A + XI 
is positive for a sufhciently large real number A, i.e all its elements are positive. 
By the Perron-Frobenius theorem |16j it has a unique eigendirection spanned 
by a positive vector. Let p be an eigenvector of this type with components 
Pi. The corresponding eigenvalue is positive. Let it be denoted by /3. Another 
consequence of the Perron-Frobenius theorem is that all other eigenvalues of 
A + XI have modulus smaller than /?. In particular the real part of any other 
eigenvalue is smaller than /3. The vector p is an eigenvector of A with eigenvalue 
a — /3 — X and all other eigenvalues of A have real part smaller than a. 

If A is a matrix of the above special form then it can be shown that the 
matrix B = A~^ is a positive matrix. One way of proving this as follows. Let 
a: be a vector in R" and consider the equation Ax = y. Inverting the matrix 
is equivalent to solving this equation for x. The equation can be written in 
components as 

auXi + tti^i-iXi-i = TJul < i < n (B.l) 
where the indices are to be interpreted modulo n. Hence 

ai+iAXi = (-ai+i,i+i)xi+i +yi+i. (B.2) 

Note that the coefficients in this equation are positive. By substituting these 
equations into each other successively with i increasing from one to n it is 
possible to obtain an equation of the form; 

^]^ai,i+i^ xi = ^]^(-aii)^ Xn +^Ciyi (B.3) 

where the coefficients Ci are positive. Rearranging this gives 

(-l)"+i(detA)x„=^c,y,. (B.4) 

Any other Xi can be determined in an analogous way. The determinant of A is 
equal to H,; + (—1)"^^ Hi Oi.j+i- Putting these facts together gives the proof 
of the desired statement. 



36 



It can be concluded from the Perron-Frobenius theorem that B has a unique 
eigendirection spanned by a positive vector. This is also an eigenvector of A 
and so must be proportional to p. The corresponding eigenvalue is and is 
positive. Hence a is positive. 
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